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SECTION  1 
INTRODUCTION 

The  finite  element  method  is  a  very  simple  way  to  find  many  things  including 
stresses  and  deflections  in  a  structure.  The  meshing  choice,  element  formulation 
and  density  determine  the  results  for  a  given  structure.  Additionally, 
computational  costs  are  very  sensitive  to  the  degrees  of  freedom  in  the  mesh. 
Beyond  those  choices,  the  integration  scheme  and  the  integration  network  also 
have  a  large  influence  in  the  results.  These  sensitivities  are  exaggerated  in 
certain  load  cases,  especially  in  a  shear  or  bending  dominated  loading 
conditions. 

This  report  is  a  simple  introduction  to  these  factors  and  the  mathematical  theory 
behind  them.  The  theory  behind  the  report  is  that  an  analyst  who  understands 
when,  as  well  as  why  each  problems  occurs,  will  be  able  to  judge  when  it  is  best 
to  select  different  options  in  the  analysis. 

1.1  Background 

The  finite  element  method  is  based  on  using  the  constitutive  laws  and  stress- 
strain  relationships  for  a  given  object  to  approximate  the  stress  field  throughout 
the  object.  Depending  on  the  element  this  is  done  in  different  ways.  Each 
element  is  given  different  degrees  of  freedom  that  affect  how  the  energy  of  the 
material  is  represented. 

However  the  energy  is  represented,  the  solution  method  is  similar.  The  potential 
energy  is  calculated  by  subtracting  the  work  done  by  the  applied  loads  from  the 
strain  energy  of  the  material.  Setting  the  variation  in  the  potential  energy  to  zero 
defines  the  equilibrium  position.  By  creating  as  many  equations  as  there  are 
degrees  of  freedom,  the  value  of  each  degree  of  freedom  can  be  found  by 
solving  the  simultaneous  equations.  The  general  principle  is  shown  below. 


k=V-W 
K =u  Ku—u  F 
8k  =  8uT Ku-8uT  F 
8k  =  0  =  SuT  \Ku  -F] 

Ku  =  F 
u  =  K~lF 

The  key  to  this  is  that  the  work  must  be  expressed  as  forces  applied  at  discrete 
points  on  the  structure  and  that  the  strain  energy  can  be  expressed  as  some 
stiffness  relationship  multiplied  by  a  combination  of  two  displacements  at  the 
nodes.  Fortunately,  the  constitutive  laws  and  strain  displacement  relationships 
most  commonly  used  apply  exactly  to  this  form. 


k  -  Total  potential  energy 

V  -  Strain  energy 

W-Work 

K  -  Stiffness  matrix 

F  -  Force  vector 

u  -  Nodal  displacements 

8-  First  Variation 
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There  are  some  flaws  with  the  method.  The  constitutive  laws,  in  certain  cases, 
create  a  situation  where  the  solution  to  the  system  of  equations  does  not 
converge  to  the  correct  answer.  There  are  a  variety  of  simple  solutions  to  rectify 
the  problem,  either  based  on  changing  the  properties  of  the  element,  or  by 
changing  the  way  the  stiffness  matrix  is  formulated. 

1.2  Scope 

This  report  will  use  as  an  example  a  cantilevered  beam  with  a  tip  shear  load 
shown  in  Figure  1 .  This  condition  was  chosen  because  it  is  relatively  easy  to  find 
a  closed  form  solution  to  the  deflection  and  that  the  problems  in  the  elements  will 
be  easily  identified.  Additionally,  the  cantilevered  beam  is  an  elementary  model 
for  many  objects  in  structural  analysis,  from  wings  to  buildings. 


Figure  1.  Test  Configuration. 

The  elements  that  are  examined  are  shell  and  brick  elements,  both  with  nodes 
only  at  the  corners.  During  the  discussions  the  equations  will  be  developed  for  a 
two-dimensional  version  of  the  element  instead  of  the  three-dimensional  ones 
used  in  the  actual  analysis.  This  is  done  only  to  significantly  reduce  the  size  of 
the  problems.  The  mathematics  is  the  same  for  the  three-dimensional  case,  but 
the  three-dimensional  case  would  have  over  twice  the  degrees  of  freedom. 

This  report  is  intended  only  as  an  introduction  to  the  mathematical  reasoning 
behind  the  finite  element  method.  It  is  intended  only  to  allow  the  user  to  decide 
when  it  is  justified  to  use  a  given  mesh  density,  a  particular  element,  or 
integration  technique.  Hopefully  it  will  allow  an  analyst  to  better  balance  the 
benefits  of  a  highly  refined  mesh  and  the  necessity  to  reduce  computational 
costs. 
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SECTION  2 
EXACT  SOLUTION 


This  particular  test  configuration  was  chosen  because  it  not  only  produces  the 
phenomenon  that  we  want  to  explore,  but  because  it  is  relatively  easy  to  find  the 
exact  form  of  the  solution.  This  will  prepare  a  baseline  to  which  the  finite  element 
models  can  be  compared. 

Two  different  solutions  are  examined.  The  first  is  the  deflection  solution,  which 
accounts  for  both  the  bending  deformation  and  the  shear  deformation.  This 
method  uses  the  variation  of  total  potential  energy  to  determine  the  governing 
equations.  The  second  is  the  frequency  analysis,  which  uses  Euler-Bernoulli 
assumptions.  These  assumptions  should  be  valid  for  the  natural  frequencies 
because  natural  frequencies  are  independent  of  loading  conditions,  so  the  shear 
loading  should  not  matter  as  much  in  the  deflection  case.  Further,  the  beam  is  a 
long  slender  beam  so  the  Euler-Bernoulli  assumptions  should  be  valid  under  any 
loading  condition.  This  method  uses  the  summation  of  forces  to  determine  an 
equilibrium  condition. 

2.1  Deformation  Solution 

The  deflection  solution  includes  both  the  bending  and  the  shear  deformation 
energy.  The  beam  is  isotropic  and  constant  cross-section.  The  dimensions  of 
the  beam  are  shown  in  Figure  2  below.  The  curvature  is  defined  as  the  change  in 
rotation  and  the  shear  is  defined  as  the  difference  between  the  slope  and  the 
rotation.  The  total  potential  energy  is  the  bending  energy  and  the  shear  energy 
less  the  work  done  by  the  tip  load. 


Figure  2.  Dimensions  of  Test  Beam. 

There  are  several  material  properties  that  must  also  be  known: 


E  -  Young’s  modulus 
G  -  Shear  modulus 
v  -  Poisson’s  ratio 
(Kx)  -  Rotation  field 
u(x)  -  Displacement  field 


k-  Shear  stiffness 
I  -  Moment  of  inertia 
P  -  Applied  tip  load 
7i  -  Total  potential  energy 
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The  first  step  is  to  simplify  the  equation  by  making  it  non-dimensional. 


1  f  vt(  d<t>  t  ,  lf^  i 

7T  =  —  £7  —  dk  +  —  Gfd - & 

2  J  I  A  9  J  I  Wv 
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l  <Lc 
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Then,  the  first  variation  is  set  to  zero  in  order  to  find  the  equilibrium  conditions. 
Integration  by  parts  allows  grouping  on  the  variation  of  each  field. 


°=J 

OL 
1 

0  =  J  <50 


f <5f  +  klj-Slj--kjS(p-k(pSlj-  +  k(j)8(j) 
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To  set  the  variation  to  zero,  each  differential  equation  multiplied  by  an  arbitrary 
variation  must  be  equal  to  zero.  The  boundary  terms  multiplied  by  an  arbitrary 
variation  must  equal  zero  at  each  boundary  as  well.  This  produces  two 
differential  equations  and  four  boundary  conditions. 


DEs : 


// 


BCs : 


k 


--W-P 


?j=i 


f  =  kQ-kj 

Co=° 


We  find  the  equilibrium  positions  for  the  unknown  displacement  and  rotation  by 
solving  the  simultaneous  differential  equations  subjected  to  the  boundary 
conditions.  This  produces  the  two  equations  for  the  solution  shown  below.  The 
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equations  are  shown  in  non-dimensional  form,  but  it  would  be  easy  to  convert 
them  to  dimensional  form. 


$  =  ~-PT+PV 

U  1  ,  1  2  P 

=  -  prf  +  w+JLri 
L  6  2  k 


2.2  Natural  Frequency 

The  natural  frequencies  are  found  by  summing  the  forces  and  moments  about  a 
typical  section  of  the  beam  shown  in  the  figure  below.  The  moment  (M)  and  the 
shear  force  (V)  cause  stresses  on  the  end  of  the  section.  We  will  employ  the 
Euler-Bernoulli  set  of  assumptions,  which  equate  the  slope  and  the  rotation. 
Because  of  that,  the  deflection  is  the  only  field  that  we  are  attempting  to  solve. 
This  is  a  dynamic  problem,  so  our  deflection,  u,  is  a  function  of  position  (x)  and  of 
time  (t). 


Figure  3.  Typical  Bending  Beam  Section. 

First  let  us  examine  the  moment  equation.  The  Euler-Bernoulli  assumptions  also 
say  that  the  rotational  inertia  is  small,  allowing  us  to  equate  the  sum  of  the 
moments  to  zero.  The  equation  is  a  first  order  approximation  of  the  system.  This 
means  that  the  moment  and  the  shear  field  can  be  looked  at  as  a  linear  function 
across  Ax.  It  also  implies  that  higher  order  terms  of  Ax  are  approximately  zero. 


A 


Vm  =M(x,t)+  (A  - / )  Ar  _  [/  (x,  t  )Ar 
ox 

dM(x,t) 
dx 

\ 

V(x,t)  = 


- V{x,t ) 


lAv  =  o 


dM  (x,  t ) 
dx 


J 


dV(x,t) 

dx 


Ax "  -  M  (x,t)  =  0 
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Next  we  sum  the  forces  using  the  same  first  order  approximations.  This  time 
equating  the  sum  of  the  forces  to  the  mass  times  the  acceleration  of  the  section. 


OX 


pAAx 


dt2 


dV(x,t ) 


Ax  =  -pAAx 


d2i 


dx  '  dt2 

We  combine  the  equations  by  differentiating  the  result  of  the  moment  equation. 


dV(x,t )  _  d2M(x,t ) 
dx  dx2 


dM(x,t) 

dx2 


The  Euler-Bernoulli  assumption  states  that  the  moment  is  equal  to  the  bending 
stiffness  (El)  times  the  second  derivative  of  the  deflection.  This  comes  from 
equating  the  rotation  with  the  slope,  meaning  that  the  curvature  is  the  second 
derivative  of  the  deflection.  We  use  this  to  reduce  our  differential  equation  to  a 
single  time  and  space  dependent  variable. 


M{x,t)=Elp \ 
dx 

d4u  pA  d2u  _ 
dx4  El  dt 2 


4  ~\2 

U  2d  u 
-  +  C  —r  =  0 


dx4  dr 

Suitable  boundary  conditions  need  to  be  developed  for  this  differential  equation. 
For  free  vibrations  we  know  that  the  free  end  will  carry  no  shear  load  or  moment. 
Using  the  earlier  definitions  for  our  shear  and  moments,  and  a  fixed  displacement 
and  slope  at  the  cantilevered  end  we  get  the  following  boundary  conditions. 


m(0,  t)  =  0 
du(o,t)  _ 
dx 


d2u(L,t) 

dx2 

d3u(L,t) 
dx 3 


Using  separation  of  variables,  with  X  being  the  separation  constant,  the  deflection 
equation  becomes: 


(p(x)  =  C[  sin(/k)+  c2  cos ((5x)+c3  sinh(/T')  +  c4  cosh  (fix) 
(54=X2c2 

After  the  boundary  conditions  are  applied,  you  reduce  the  system  to  two 
equations  and  two  unknowns. 
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sin 

cos 


(j6L)+  sinh(BL)  cos^}4^osh(/JL) 
(j6L)  + cosh (p/J  -sin(^)+sinh(j6L) 


V. 

1 

O 

1 

1 

o 

To  avoid  the  trivial  solution,  it  is  necessary  to  make  to  determinate  of  the  matrix 
equal  to  zero.  This  is  done  by  varying  the  values  of  PL.  The  first  few  values  that 
create  a  zero  determinate  are  listed  below.  From  the  time  dependent  portion  of 
the  separation  of  variables,  the  frequency  of  oscilliation  was  shown  to  be  X. 

Using  the  material  properties  and  the  relationships  shown  above  the  PL  values 
are  changed  into  frequencies.  The  units  are  then  changed  from  radians  per 
second  to  hertz. 


PL=  1.875104 
pl_=  4.694091 
pL=  7.854757 
pL=  10.995541 
pL=  14.137169 


Xi=  9.7260  Hz 
X2=  60.9518  Hz 
X3=  170.6668  Hz 
XA=  334.4388  Hz 
X5=  552.8514  Hz 


By  changing  the  orientation  of  the  beam  it  is  easy  to  find  the  lateral  bending 
frequencies.  Only  the  first  bending  is  low  enough  to  be  relavent  to  the  problem. 


7iat=  194.5200  Hz 


This  type  of  analysis  will  not  give  the  torsional  mode  frequencies.  To  get  those 
frequencies,  we  must  look  again  at  the  typical  section  of  the  beam,  seen  in 
Figure  4. 


Figure  4.  Typical  Torsional  Beam  Section. 


This  time  there  are  no  forces  and  the  rotations  are  not  considered  small,  so  the 
moments  are  set  equal  to  polar  moment  multiplied  by  the  angular  acceleration. 

3MM-p(/  +/  =  o 

at  "  "dr 

The  moment  can  also  be  written  in  the  following  form, 


2-5 


where  J  is  the  torsional  stiffness  of  the  beam.  For  a  beam  with  a  10:1  aspect 
ratio  cross-section,  such  as  the  one  in  the  model,  J  is  defined  as  below. 

As  in  the  bending  case,  this  relationship  is  used  to  create  a  single  field  differential 
equation.  The  boundary  conditions  of  a  one  fixed  end  and  one  with  no  moment 
are  also  translated  mathmatically. 

J  =  dbh3 
d  =0.313 


Using  separation  of  variables  and  the  boundary  conditions  as  in  the  bending 
case,  the  non-trivial  solution  is  satisfied  if  X  has  the  following  values. 


d2(j){x,t)  |  g2  _Q  BCs  \ 

dx2  C  dt2  (j>(0,t)=0 

2  _  Pi1  xx  +Iyy)  d(j)(L,  t )  _ 

GJ  dx 


pl_=  ti/2 
pl_=  3ti/2 
pl_=  571/2 
pL=  7n/2 
pL=  971/2 


Xl=  89.6567  Hz 
X2=  268.9700  Hz 
X3=  448.2834  Hz 
X4=  627.5968  Hz 
X5=  806.91 01  Hz 


The  first  ten  frequencies  have  been  calculated  as: 


1  bend 

2  bend 

1  tor 

3  bend 

1  lat 

9.726 

60.952 

89.657 

170.667 

194.520 

2  tor 

4  bend 

3  tor 

5  bend 

4  tor 

268.970 

334.439 

448.283 

552.851 

627.597 

To  calculate  the  mode  shapes  we  need  to  return  to  the  shape  function.  Then  we 
can  vary  pL  from  0  to  the  value  calculated  above.  We  must  also  solve  for  each 
of  the  four  constants  that  appear  in  the  equation.  Unfortunately  we  used  some  of 
the  information  to  determine  the  frequencies,  so  we  still  have  one  more  unknown 
variable  than  the  number  of  equations.  We  can  pick  a  unit  value  for  one  of  the 
constants  and  define  the  other  by  the  ratio  of  the  two.  The  value  of  the  deflection 
is  normalized  by  the  stiffness  matrix  such  that  the  mode  shape  vector  transposed 
times  the  stiffness  matrix  times  the  mode  shape  vector  is  equal  to  one.  The 
analagous  method  of  normalization  in  the  exact  form  of  the  solution  would  be  to 
integrate  the  stiffnesses  times  the  degrees  of  freedom  squared  over  the  length 
of  the  beam.  That  is  more  effort  than  the  resulting  benefit  would  yield,  so  only 
the  unnormalized  results  are  shown  in  the  figure  below.  Also  the  torsional  modes 
are  more  difficult  to  show  when  assumed  to  be  a  part  of  a  one-dimensional 
beam.  Therefore,  Figure  5  only  shows  the  first  three  bending  modes. 
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Figure  5.  Theory  Based  Bending  Mode  Shapes. 
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SECTION  3 

SHELL  ELEMENT  LOCKING 


The  first  element  that  we  will  look  at  is  the  shell  element.  For  the  mathematical 
development  of  the  element,  a  one-dimensional  model  is  used,  but  the  premise  is 
the  same  as  in  the  full  two-dimensional  case.  The  element  mode  shapes  will  be 
presented  in  the  one-dimensional  case.  The  deflection  solution  was  calculated 
by  building  the  beam  and  solving  the  analysis  in  Excel. 

3.1  Element  Mode  Shapes 

The  first  step  that  we  need  to  accomplish  when  building  the  element  is  to  define 
the  displacement  field  across  the  element.  The  displacement  at  each  node  must 
be  completely  independent  of  the  other  nodes  in  the  element.  This  is  done  by 
creating  a  shape  function  that  has  a  value  of  one  at  one  node  and  zero  at  all  the 
other  nodes.  These  shape  functions  are  added  together  in  a  linear  combination 
with  each  shape  function  multiplied  by  the  value  at  a  particular  node.  The  shape 
functions  are  written  with  respect  of  the  local  variable  r  and  can  be  seen  in  Figure 
6. 


h ,  =  —  (r  + 1) 
1  2 


u(r)  =  /?j  (r  )u  i  +  h2  ( r)u2 
(j)(r )  =  /z,  (r )0,  +  h2  {r)(j)2 


Figure  6.  Shell  Element  Shape  Functions. 

The  displacement  field  can  also  be  written  in  matrix  form.  This  will  be  a  more 
advantageous  form  for  later  work. 


W, 

u{r)  _ 

\ 

0  !  h2  0  " 

0i 

_0 

\  !  0  h2_ 

u2 

02 

u(r)  =  H(r)u 


The  next  step  is  to  define  the  strain  field  in  the  element.  There  are  two  types  of 
strain  just  as  it  was  in  the  development  of  the  exact  solution.  The  bending  strain 
and  the  shear  strain  are  defined  the  same  way  as  before. 

^_d(j)  _  cl(j)  dr  _  2  d(j) 

dx  dr  dx  l  dr 

du  ,  du  dr  ,  2  du 

/  = - 0  = - 0  = - 0 

dx  dr  dx  l  dr 
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The  curvature  and  the  shear  deformation  can  also  be  expressed  in  matrix  form. 
The  B  matrix  is  called  the  strain-interpolation  matrix.  In  this  notation,  the  prime 
marking  indicates  a  derivative  with  respect  to  the  local  coordinate  r,  and  I  is  the 
global  element  length. 


r 

2  7 

2  7 

W, 

K 

0 

Yh' 

o  -K 

l  2 

7 

2  ,  / 
Ylh' 

—  h'0  -  h7 

1  2  2J 

u2 

02 

e(r)  =  B(r)il 


Expressing  the  our  strain  energy  equation  similarly  to  the  exact  solution  form, 
and  then  converting  it  into  matrix  notation  points  out  the  form  for  the  stiffness 
matrix  described  in  the  background  portion  of  the  introduction.  J  is  the  Jacobian 
matrix  expressing  the  transformation  from  the  local  to  the  global  coordinates.  In 
this  case  there  is  only  one  transformation,  which  was  earlier  expressed  as  1/2. 


V  =E\eI(KY  Fir  +1-\gk(j) 


i  l 


dr 


-1 


1  f  1 

'El  0  " 

K 

v=yfl_K  rj 

. 0  GK 

7 

V  = 


1 

—  u 
2  “ 


\BTCB\\j\\lr 


—  dr 
2 


'*'•  f  A  T  t y 

u  =  —u  Ku 
~  2  “  “ 


The  stiffness  matrix  is  a  way  of  expressing  the  relationship  between  an  applied 
force  at  one  position  on  a  structure,  and  the  response  at  some  other  position. 
The  relationships  are  assumed  to  be  linear,  letting  us  use  the  follow  form  for  the 
equilibrium  position.  Linearity  also  allows  us  to  combine  several  elements 
together  by  adding  the  common  matrix  entries  for  common  nodes. 


Ku  =  F 


There  are  several  interesting  aspects  to  a  structure’s  stiffness  matrix.  The  mode 
shapes  that  an  element  or  structure  can  assume  and  their  frequencies  can  be 
found  by  extracting  the  eigenvalues  and  eigenvectors  from  the  stiffness  matrix. 
The  eigenvectors  describe  the  displacement  of  each  node  while  the  associated 
eigenvalue  predicts  the  natural  frequency  for  that  mode. 

The  shell  element  has  four  modes,  two  rigid  body  modes  and  two  elastic  modes. 
The  rigid  body  modes  have  a  frequency  of  zero,  while  the  elastic  modes  have 
non-zero  frequency.  The  eigenvectors  are  shown  with  normalized  values  that 
make  the  magnitude  of  the  vector  one.  It  is  important  to  also  note  that  each 
mode  is  orthogonal  to  the  others.  Any  normalized  eigenvector  multiplied  by 
another  is  equal  to  zero,  while  any  normalized  eigenvalue  squared  is  equal  to 
one. 
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Figure  7.  Shell  Element  Mode  Shapes. 


3.2  Deflection  Solution 

The  final  solution  is  a  linear  combination  of  the  mode  shapes;  each  multiplied  by 
a  modal  participation  factor.  It  is  easier  to  solve  the  set  of  simultaneous  linear 
equations  than  it  is  to  find  the  modal  participation  factors.  Most  times  the  modal 
participation  factors  are  extraneous  to  the  final  deflection  and  stress  results. 

Demonstrating  the  ease  of  using  FEM  technique  on  simple  geometry,  the  system 
is  solved  for  using  Excel.  The  speed  that  the  solution  approaches  the  exact 
solution  is  dramatic.  Unfortunately,  because  of  the  shear  locking,  it  will  not 
become  apparent  until  changes  are  made  to  the  integration  scheme. 
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The  meshing  scheme  used  on  the  beam  is  to  divide  it  equally  along  the  length  of 
the  beam.  Cases  of  one,  two,  and  three  elements  are  inspected. 


Gaussian  integration  is  used  to  evaluate  the  stiffness  matrix.  By  sampling  a 
limited  number  of  places,  you  can  evaluate  a  polynomial  integral  exactly  the  over 
a  region.  The  number  of  sampling  points  is  determined  by  the  degree  of  the 
polynomial  that  you  wish  to  evaluate.  N  sample  points  can  exactly  evaluate  an 
integral  of  (2n  -1 )  order.  The  polynomial  in  the  stiffness  matrix  is  second  order 
so  there  must  be  a  minimum  of  two  sample  points.  The  points  must  be  the 
zeroes  of  the  Legendre  Polynomial.  For  n  equal  two  this  means  that  the  sample 
points  will  be  r=±3'&.  There  must  also  be  a  weighting  factor  multiplied  by  the 
results  at  each  sample  points.  For  n  equal  two,  the  weight  factor  happens  to  be 
one  for  both  points.  This  yields  the  following  as  the  integrated  stiffness  matrix  for 
an  element  length  of  ten. 


K  = 


183081 

915404 

-183081 

915404 


915404  -183081  915404 
6103177  -915404  3050863 
-915404  183081  -915404 
3050863  -915404  6103177 


The  first  two  degrees  of  freedom  are  constrained  to  zero  displacement.  This 
means  when  solving  for  the  equilibrium  solution  it  is  only  necessary  to  solve  for 
the  last  two  degrees  of  freedom.  The  applied  load  is  a  point  force  on  the 
unconstrained  node.  That  makes  our  final  system  of  equations  the  following. 
The  same  is  done  for  the  two  and  three  element  systems.  The  difference  being 


u2 

'  183081 

-915404" 

-1 

100 

$2 

-915404 

6103177 

0 

the  common  degrees  of  freedom  are  added  together  creating  a  banded  stiffness 
matrix.  The  displacements  are  shown  in  the  following  chart  compared  to  the 
exact  displacements.  _ _ 


Number  of 
Elements 

Tip  Displ. 

Exact 

6.897098 

One 

Element 

0.002184 

Two 

Elements 

0.008729 

Three 

Elements 

0.019609 

Figure  8.  Locked  Beam  Displacements. 
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It  is  obvious  from  either  Figure  7,  or  the  table  of  tip  displacements  that  the  locking 
phenomenon  is  a  severe  flaw  in  the  method.  It  is  clear  something  is  very  wrong 
in  this  simple  case,  but  in  complex  loading  of  a  shell  modeled  blade  it  would  not 
be  as  easy  to  determine  that  the  shear  deformation  is  not  contributing. 

3.3  Cause  of  Locking 

An  over-constrained  system,  or  perhaps  better  put,  an  overly  inflexible  set  of 
shape  functions,  causes  shears  locking.  Two  linear  functions  only  have  four 
unknowns  to  determine.  After  applying  the  four  boundary  conditions  there  are  no 
more  unknowns  to  manipulate.  This  means  that  the  chances  that  the  shape 
functions  will  satisfy  the  differential  equations  are  slim.  It  is  easy  to  see  this 
through  working  the  example  we  have  been  working. 

~(7?)  =  ci77  +  c2  (j)(ri)=c3ri  +  c4 

u  ,  X  0,(t)  =  c3 

TM=c' 

The  finite  element  code  does  not  apply  all  the  boundary  conditions  and  then  try 
to  satisfy  the  differential  equations.  It  applies  the  geometric  boundary  conditions, 
then  lets  the  natural  boundary  conditions  come  through  the  minimization  process. 
Following  this  procedure,  we  start  by  applying  the  geometric  boundary 
conditions. 


*1,.,  =»=<■, 


77=0 


)=c1ri 


Next,  the  differential  equations  are  satisfied.  This  produces  the  locking  results. 


c3  -  0  =  0  =>  c3  =  0 

0(77)= 0 


/ 

k$-f-k^-  =  0 

0 — 0 — kc ,  =  0  Cj  =0 


u 

L 


(77)  =  0 


The  actual  results  are  not  zero  displacement  because  the  minimization  process 
tries  to  impose  the  natural  boundary  conditions  while  solving  the  differential 
equations.  This  would  enforce  some  tip  displacement,  but  it  can  be  shown  that 
the  tip  displacement  approaches  zero  as  the  aspect  ratio  of  the  element 
increases. 


3.4  Element  Applicability 

Most  commercial  packages  do  not  use  the  type  of  integration  scheme  on  shell 
elements  because  of  the  severe  problems  shown  above.  One  notable  exception 
is  when  analyzing  pressure  vessels.  Some  codes  use  the  fact  that  the  fully 
integrated  shell  element  does  not  react  to  bending  loads  to  better  model  the 
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steady  state  equilibrium  of  pressure  vessels.  The  elements  will  only  deform 
axially,  so  surfaces  will  remain  perpendicular  to  the  pressure. 

Otherwise  codes  will  automatically  use  reduced  integration.  It  will  be  shown  later 
that  reduced  integration  is  an  easy  solution  to  the  shear-locking  problem. 
Reduced  integration  causes  some  other  problems  of  which  the  user  must  be 
aware. 

The  constitutive  law  and  strain  displace  relationships  will  enforce  equilibrium 
within  the  element  even  if  flawed  boundary  conditions  provide  the  incorrect 
forces.  The  finite  element  method  does  not  enforce  equilibrium  between  the 
elements.  That  is  why  an  increasing  number  of  elements  creates  a  better 
solution.  It  should  be  noted  that  it  is  not  the  correctness  of  the  FEM  that  causes 
a  finer  mesh  to  be  more  exact  in  this  case,  but  the  flaw  in  the  method  that  does. 


3-6 


SECTION  4 

BRICK  ELEMENT  LOCKING 


Like  the  shell  element  development,  the  brick  element  will  also  be  pursued  one 
dimension  lower  than  the  actual  element.  The  brick  element  will  be  shown  as  a 
two  dimensional  element  like  a  membrane  element.  This  reduces  the  elemental 
stiffness  matrix  size  from  24x24  to  8x8.  The  modes  shapes  are  calculated  on 
Maple,  while  the  deflection  solution  was  found  using  Nastran. 

4.1  Potential  Energy 

The  potential  energy  equation  is  very  different  from  the  ones  used  before.  This  is 
because  the  dimension  of  depth  is  added,  eliminating  the  need  to  add  redundant 
degrees  of  freedom  for  rotation.  The  potential  energy  is  defined  by  the  Cauchy 
formulation  of  strain  energy.  Strain  energy  is  likened  to  work  done  on  the 
material  that  will  be  recovered  when  the  load  is  removed. 


Figure  9.  Brick  Element  Dimensions 

4.2  Element  Mode  Shapes 

Again,  the  mode  shapes  start  with  the  shape  functions  for  the  element.  The  local 
coordinates  used  to  define  the  shape  functions  extend  from  negative  one  to  one 
in  each  direction.  The  nodal  numbering  is  defined  starting  from  the  first  quadrant 
increasing  by  the  right-hand  rule. 


s 


-1  1 


Figure  10.  Brick  Element  Local  Coordinates. 
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The  shape  functions  have  a  unit  value  at  the  node  of  interest,  and  zero  values  at 
all  other  nodes.  A  typical  shape  function  is  shown  after  the  shape  functions 
below. 


K  =^(l  +  r)(l  +  s) 


h3  =^(1-?')(l-^) 

K  =^(|  +  r)(|-0 


Figure  11.  Brick  Element  Shape  Function. 


The  next  step  is  to  define  the  strain-displacement  relationship.  The  differences 
from  the  shell  element  are  axial  strain  in  two  dimensions  and  that  the  shear  strain 
is  defined  differently.  The  k  index  in  the  matrix  notation  is  to  denote  the  shape 
function  number.  The  strain  interpolation  matrix  will  be  a  three  by  eight  matrix  in 
this  example. 
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The  final  ingredients  needed  to  find  the  stiffness  matrix,  and  then  the  mode 
shapes,  are  the  constitutive  laws.  The  plane  stress  relationships  are  used  to 
define  the  connection  between  the  strain  and  the  stress. 
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The  Jacobian  becomes  a  more  complex  matrix  because  of  the  second  dimension 
added  to  the  problem.  The  Jacobian  is  defined  in  the  following  manner.  The 
Jacobian  is  used  to  transform  the  coordinates  from  the  global  to  the  elemental 
coordinates. 


dx 

dx 

dx 

dr 
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dx 
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Combining  all  the  relationships  into  the  strain  energy  equation  will  produce  the 
stiffness  matrix,  just  as  it  did  for  the  shell  element.  The  same  notation  applies  as 
before.  The  B  matrix  is  the  strain  interpolation  matrix.  The  matrix  used  in  the 
constitutive  law  is  the  C  matrix. 
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The  stiffness  matrix  for  this  element  is  an  eight  by  eight  matrix.  This  means  there 
will  be  eight  mode  shapes.  There  are  three  rigid  body  modes,  two  translations 
and  one  rotation.  There  are  five  elastic  modes,  one  uniform,  two  shear,  and  two 
bending  modes.  All  of  the  modes  are  shown  below  with  the  associated 
eigenvalue.  The  rigid  body  modes  have  a  zero  eigenvalue,  meaning  that  they  all 
must  be  constrained  or  the  stiffness  matrix  will  be  singular. 


Horizontal  Translation 

X=0 


Vertical  Translation 

X=0 


Rotation 

X=0 


Figure  12.  Brick  Element  Rigid  Body  Modes. 
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Bending 

^=0.4976 


Bending 

^=0.4976 


Uniform  Contraction 
?i=1 .47 


X-0.7575 


Figure  13.  Brick  Element  Bending  Modes. 


4.3  Deflection  Solution 

Patran  generated  the  code  to  analyze  the  deflection  solution  in  Nastran.  There 
were  two  models  used.  The  first  model  had  180  elements,  each  with  eight 
nodes.  The  elements  each  had  a  length  and  width  of  a  third  and  were  one-tenth 
inch  deep.  This  model  is  shown  in  Figure  14.  The  other  model  had  the  same 
dimensions,  but  the  elements  were  one-fifth  inch  long  and  wide. 


Figure  14.  Brick  Element  Beam  Model. 

Both  beams  are  held  to  the  same  boundary  conditions.  At  one  end  of  the  beam 
the  end  nodes  were  set  to  zero  displacement  out  of  the  plane.  The  center  two 
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nodes  were  additionally  held  to  no  displacement  within  the  plane  either.  The 
other  end  had  a  100-pound  downward  force  applied  to  the  upper  center  node. 

There  are  several  different  integration  methods  that  can  be  used  by  Nastran  to 
evaluate  the  stiffness  matrix.  The  option  used  in  this  model  told  Nastran  to  use 
full  integration  to  evaluate  the  stiffness  matrix.  This  meant  for  two  gauss  points 
to  be  used  in  each  direction  in  obtaining  the  exact  solution  since  one  linear 
equation  multiplied  by  another  creates  a  quadratic  equation.  As  in  the  shell 
element,  full  integration  under  shear  loading  creates  locking. 


Figure  15.  Locked  Brick  Element  Deflection  Solution. 

Even  with  fifty  elements  along  the  span  of  the  beam,  the  deflection  is  only  about 
one-third  of  what  the  exact  solution  is.  This  shows  that  although  the  model  will 
still  converge  to  the  correct  solution  through  the  discontinuity  of  equilibrium,  the 
efficiency  of  the  method  is  dramatically  reduced. 

4.4  Natural  Frequencies 

Using  the  same  integration  scheme,  but  using  a  modal  analysis  instead  of  a 
linear  analysis,  Nastran  will  give  the  first  ten  modes  and  their  associated 
frequencies.  These  are  compared  to  the  calculated  values  in  the  table  below. 


Calculations  Error 
-3.11% 


1  lateral 

1  torsion 

2  torsion 

3  torsion 

4  torsion 

1  bend 

2  bend 

3  bend 

4  bend 

5  bend 


FEM 

188.48 

99.24 
311.76 

563.24 
872.36 
16.34 
102.29 
285.75 
559.20 
922.83 


194.52 

89.66 

268.97 

448.28 

627.60 

9.73 

60.95 

170.67 

334.44 

552.85 


10.69% 

15.91% 

25.64% 

39.00% 

67.97% 

67.82% 

67.43% 

67.21% 

66.92% 


The  locking  effect’s  dependence  on  the  aspect  ratio  of  the  element  and  the 
structure  is  clear.  In  the  lateral  deflection,  the  structural  length  (10)  divided  depth 
(2)  is  five,  while  in  the  bending  direction,  the  quotient  (10/0.1)  is  100.  The  error  is 
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approximately  twenty  times  greater  as  is  the  quotient.  Further,  the  element  has 
an  aspect  ratio  of  over  three  to  one  for  the  bending  mode  and  just  one  to  one  for 
the  lateral  bending. 

The  percent  error  is  also  dependent  on  the  mode  shape  number.  Some  modes 
have  different  amounts  of  curvature  in  each  direction.  The  more  curvature  that  is 
required,  the  more  error  there  will  be  in  general.  The  bending  modes  are  all 
similar  in  error,  potentially  because  it  is  the  maximum  error  limit. 

4.5  Cause  of  Locking 

The  element  locks  because  it  is  an  inflexible  system  just  as  the  shell  element  is. 
The  definition  of  linear  is  a  little  different  in  these  shape  functions,  creating  a 
more  flexible  system,  but  even  this  added  flexibility  is  not  enough  to  unlock  the 
system.  The  existence  of  four  nodes  implies  that  we  will  need  four  coefficients  to 
correctly  determine  the  system.  Unfortunately,  two  linear  functions  of 
independent  variables  only  have  three  coefficients  when  combined  since  the 
constants  can  be  added  together.  This  is  how  a  non-linear  term  gets  added  into 
the  linear  shape  functions  shown  below. 

u(x,  y)  =  clx  +  c2y  +  c3xy  +  c4 
v(x,  y)  =  c5x  +  c6y  +  c1xy  +  ci 

Just  as  in  the  example  with  the  shell  element  the  brick  element  is  subject  to  both 
geometric  and  natural  boundary  conditions.  Each  displacement  field  has  one  of 
each  type  at  both  ends  of  the  element.  Take  as  an  example  the  element  with  the 
coordinate  axes  shifted  to  the  position  shown  below  to  make  the  equations  more 
easily  examined.  The  single  equation  requiring  no  displacement  in  the  u 
direction  acts  as  two  separate  constraints  at  x=0. 

w(0,  y)  =  c2y  +  c4  =0 
v(0,y)  =  c6y  +  c8  =0 

but  because  the  equation  must  hold  for 
all  values  of  y  the  only  solution  is: 

c2  =  0  c4  =  0 

c6=0  c8  =0 

Figure  16.  Brick  Element  with  Shifted  Axes. 

Likewise,  the  boundary  conditions  on  the  other  end  relating  deflection  to  the 
applied  shear  also  add  four  constraints.  These  eight  constraints  combined  with 
two  differential  equations  create  a  total  of  ten  constraints  and  only  eight 
unknowns.  The  system  is  over-constrained  and  tends  toward  the  trivial  solution 
trying  to  satisfy  all  conditions. 
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SECTION  5 

SHELL  ELEMENTS  WITH  REDUCED  INTEGRATION 

Locking  needs  to  be  removed  in  order  to  predict  the  stress  correctly  and 
displacement  fields  in  some  structures.  There  must  be  changes  made  to  the 
element  to  have  it  work  correctly.  There  are  two  main  ways  to  this:  change  the 
degrees  of  freedom  to  match  the  number  of  constraints  or  modify  how  the 
stiffness  matrix  is  built  so  that  the  constraints  can  be  followed  with  the  same 
number  of  degrees  of  freedom. 

The  number  of  degrees  of  freedom  is  dependent  on  the  number  of  nodes.  By 
introducing  a  node  at  the  mid-span  point  of  the  element,  two  additional  degrees 
of  freedom  are  introduced.  This  brings  the  total  to  three  translations  and  three 
rotations,  equaling  the  constraints.  There  is  no  locking  in  this  new  element.  This 
is  similar  to  the  method  followed  in  the  brick  element  with  a  bubble  mode. 
Changing  the  shape  functions  to  a  higher  order  polynomial  increases  the  number 
of  constants  that  must  be  solved  for  to  at  least  equal  to  the  number  of 
constraints. 

We  will  leave  that  for  the  next  section  and  now  see  what  can  be  done  if  there  is  a 
need  to  use  a  two-node  element  instead  of  the  one  with  three  nodes.  If  the 
stiffness  matrix  is  determined  by  using  a  lower  order  Gaussian  integration  the 
element  will  no  longer  have  locking.  This  section  will  focus  on  this  reduced  order 
integration  technique. 

5.1  Mode  Shapes 

The  shape  functions  will  still  be  linear  for  this  analysis.  This  makes  the  function 
to  be  integrated  a  second  order  polynomial  when  the  shape  functions  are 
multiplied  together.  To  exactly  evaluate  a  second  order  polynomial  would  require 
two  Gauss  points.  For  this  analysis  we  only  use  one  Gauss  point. 

This  changes  the  values  of  the  stiffness  matrix  and  therefore  possibly  the 
eigenvalues  and  eigenvectors.  The  overall  eigenvectors  do  not  change,  but  an 
eigenvalue  does  change.  The  bending  mode  has  an  eigenvalue  that  was 
lowered  by  three  orders  of  magnitudes.  This  means  that  the  amount  of  energy 
that  is  needed  to  cause  a  change  in  the  deflection  is  much  less  than  the  locking 
situation,  allowing  the  locked  configuration  to  unlock. 

The  reduced  order  integration  will  exactly  integrate  the  equations  in  the  stiffness 
matrix,  except  for  those  in  the  second  or  fourth  column  and  the  second  or  fourth 
row.  That  is  because  when  the  strain  interpolation  matrix  is  transposed  and 
multiplied  by  itself,  only  those  positions  have  a  term  of  r  squared.  This  is 
apparent  when  we  recall  the  strain  interpolation  matrix  that  was  developed  before 
is  as  shown  below. 
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The  only  mode  shape  that  is  completely  defined  by  the  stiffness  in  the  rotations  is 
the  fourth  mode.  It  is  also  the  only  mode  that  changes  eigenvalues.  All  other 
eigenvalues  are  based  on  the  stiffness  of  either,  combined  rotation  and 
translation  (odd  column  and  even  row  and  vice  versa),  or  pure  translation  (odd 
rows  and  columns). 

5.2  Force-Displacement  Matrix 

The  construction  of  the  stiffness  matrix  through  full  integration  of  the  constitutive 
equations  and  the  strain  displacement  relationship  effectively  minimizes  the 
unknown  constants  in  the  deflection  solution  in  a  given  ratio  in  accordance  with 
the  governing  laws.  The  matrix  inversion  and  the  applied  forces  then  solve  for 
the  magnitude  of  that  ratio  by  creating  a  linear  combination  of  those  deflection 
ratios.  This  is  obvious  when  you  view  the  equation. 

d  =  k~'f 


0.001912  0.001147 
0.001147  0.000689 
0.005732  0.003441 
0.001147  0.000689 
0.009552  0.005736 
0.001147  0.000689 


0.005732  0.001147  0.009552 
0.003441  0.000689  0.005736 
0.019103  0.004589  0.034383 
0.004589  0.001378  0.009177 
0.034383  0.009177  0.066854 
0.004589  0.001378  0.010324 


0.001147 

0.000689 

0.004589 

0.001378 

0.010324 

0.002067 


The  inverse  of  the  stiffness  matrix  is  clearly  a  way  to  combine  the  individual 
shape  functions.  To  represent  this,  the  figure  below  shows  the  deflection  for  a 
unit  force  applied  at  a  particular  node.  The  final  solution  is  a  linear  combination 
of  each  of  these  shapes,  each  with  the  participation  coefficient  equal  to  the  load 
applied  at  the  respective  node. 


Figure  17.  Beam  Deflection  Mode  Shapes. 
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The  first  thing  that  can  be  seen  is  that  the  unit  force  creates  a  greater 
displacement  when  applied  at  a  more  outboard  node.  This  makes  sense 
because  it  would  create  a  greater  moment  about  the  beam  root. 

The  inverted  stiffness  matrix  also  shows  that  the  strain-displacement 
relationships  are  enforced.  This  is  best  seen  in  the  first  column  representing  a 
load  applied  at  the  first  node.  The  rotation  at  all  further  nodes  are  equal  to  the 
rotation  at  the  first  node.  This  is  because  the  beam’s  curvature  is  proportional  to 
the  moment  at  that  point.  Outboard  of  the  applied  force  there  is  no  reaction 
force,  nor  reaction  moment.  This  means  that  curvature  outboard  should  be  zero 
implying  a  constant  rotation  and  slope.  Both  of  these  implications  are  found  to 
be  true  in  the  graph  above. 

5.3  Deflection  Solution 

The  difference  that  this  makes  in  the  final  solution  is  drastic.  The  effect  of 
making  the  change  in  rotation  easier  makes  the  solution  converge  towards  the 
exact  solution  very  quickly. 

Once  again  the  solution  was  done  on  Excel.  The  stiffness  matrix  was  inverted  by 
using  the  row-echelon  format.  The  inverted  matrix  multiplied  by  the  force  vector 
will  produce  the  deflection  solution.  The  solution  is  found  for  one,  two,  and  three 
elements  across  the  length  of  the  beam.  Solutions  are  compared  both  in  the 
graph  and  in  the  table  of  tip  displacements. 


Number  of 
Elements 

Tip  Disp. 

Exact 

6.897098 

One 

Element 

5.17296 

Two 

Element 

6.466063 

Three 

Element 

6.705527 

Figure  18.  Reduced  Integration  Shell 
Deflection  Solution. 


The  elements  are  more  accurate  as  the  actual  deflection  approaches  a  linear 
function.  When  there  is  more  curvature  in  the  solution  the  elements  must  be 
smaller  to  stay  close  to  the  exact  answer.  Even  so,  the  error  is  small  considering 
there  are  only  three  elements.  In  fact,  the  solution  converges  to  the  exact 
solution  very  quickly.  Remember  that  the  locked  solution  yields  a  zero  deflection 
solution  for  high  aspect  ratio  elements  such  as  these. 
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The  effect  the  change  has  on  the  rotation  is  more  curious.  The  reduced  energy 
to  cause  rotation  allows  the  tip  rotation  to  equal  the  exact  for  any  number  of 
elements,  including  only  one.  Therefore  the  stiffness  matrix  applies  the  tip 
moment  boundary  condition  exactly. 

5.4  Natural  Frequencies 

To  find  the  first  ten  natural  frequencies  of  the  system  there  must  be  at  least  ten 
degrees  of  freedom.  This  means  that  you  would  need  to  invert  at  least  a  ten  by 
ten  matrix,  more  if  you  wish  to  have  more  accuracy.  This  makes  the  reverse 
echelon  method  impractical.  Instead,  Nastran  solved  the  system  of  equations. 

For  this  example  a  system  of  80  four-node  shell  elements  predicted  the 
frequencies.  The  beam  had  four  elements  across  the  width  and  twenty  along  the 
length.  A  figure  showing  the  mesh  is  located  below. 


Figure  19.  Meshing  of  Reduced  Order  Shell  Elements. 


The  higher  number  of  elements  attempts  to  keep  the  curvature  close  to  zero 
across  any  given  element.  This  will  reduce  the  error  in  energy  in  a  given  mode 
and  make  the  frequency  closer  to  the  exact  frequency.  The  improvement  of  the 
frequency  prediction  over  the  locked  brick  elements  is  apparent.  These  could  be 
improved  by  adding  more  elements.  They  would  be  best  added  across  the  width 
to  rectify  the  error  in  the  torsion  modes. 


FEM 

Calculations 

Error 

1  lateral 

188.42 

194.52 

-3.14% 

1  torsion 

90.79 

89.66 

1 .26% 

2  torsion 

276.92 

268.97 

2.96% 

3  torsion 

476.23 

448.28 

6.23% 

4  torsion 

696.26 

627.60 

10.94% 

1  bend 

9.85 

9.73 

1 .23% 

2  bend 

61.50 

60.95 

0.90% 

3  bend 

172.12 

170.67 

0.85% 

4  bend 

337.46 

334.44 

0.90% 

5  bend 

557.91 

552.85 

0.98% 

The  higher  order  torsion  modes  have  the  most  error  because  they  have 
components  of  curvature  in  both  the  length  and  width  directions.  The  higher  the 
mode  the  more  nonlinear  the  deflection  is  across  a  given  element,  causing  the 
error  to  increase  as  the  modal  number  increases. 
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5.5  Element  Applicability  and  Meshing  Criteria 

This  element  is  a  good  one  to  use  in  most  situations.  It  is  the  only  integration 
technique  used  by  most  codes  for  shell  elements.  The  major  limitation  to  the 
reduced  integration  techniques  is  that  it  cannot  be  used  when  the  stiffness  or 
mass  matrix  is  not  symmetric.  So,  if  there  is  an  active  control  system  or  energy 
losses  not  associated  with  damping  forces  the  reduced  integration  technique  will 
not  work. 

The  best  method  of  meshing  is  to  have  as  many  elements  as  possible  across  the 
stress  gradient,  the  steeper  the  stress  gradient  the  more  elements  required.  In 
our  example  the  stress  gradient  lines  run  across  the  width  of  the  beam.  This 
implies  that  for  best  results  single  elements  could  run  across  the  width  of  the 
beam,  giving  the  most  possible  elements  across  the  stress  flow. 

In  most  cases  the  stress  gradient  is  not  known  before  the  analysis.  The  best 
idea  is  to  then  keep  the  elements  close  to  an  aspect  ratio  of  one  to  get  an  initial 
mesh.  After  the  general  stress  field  is  known  a  refined  mesh  can  be  created. 

There  will  always  be  a  tradeoff  between  the  increase  in  accuracy  and  the 
increase  in  computational  costs.  The  increase  of  computational  time  is  rough 
proportional  to  the  number  of  elements  squared.  Ultimately  the  solution  to  this 
tradeoff  is  left  to  the  analyst. 
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SECTION  6 

BRICK  ELEMENT  WITH  BUBBLE  NODE 

The  same  concepts  that  could  be  applied  to  the  shell  element  to  relieve  the 
locking  can  be  applied  to  the  brick  element.  Either  the  element  can  be  given 
more  degrees  of  freedom  or  the  stiffness  matrix  can  be  evaluated  a  different  way. 

This  section  will  look  at  how  the  addition  of  an  additional  node,  not  on  the 
boundary  with  another  element,  can  add  the  needed  flexibility  to  unlock  the 
system.  The  node  will  be  placed  at  the  center  of  the  element  in  the  two 
dimensional  model  of  the  full  brick  element.  The  additional  node,  with  its 
associated  degrees  of  freedom,  will  bring  the  total  degrees  of  freedom  equal  to 
the  total  number  of  constraints. 

6.1  Shape  Functions  and  Mode  Shapes 

There  will  be  changes  to  the  mode  shapes  because  of  the  introduction  of  the 
additional  node.  This  new  shape  function  should  have  a  value  of  one  at  the 
origin  and  zero  values  on  the  edges.  Additionally,  the  original  shape  functions 
must  be  changed  as  well.  They  all  had  a  value  of  one  quarter  at  the  origin, 
whereas  they  should  now  have  a  value  of  zero.  This  is  accomplished  by 
subtracting  a  quarter  of  the  new  shape  function  from  each  of  the  old  ones.  The 
shape  functions  are  shown  here. 

K  =^(l+/'Xl  +  ^)-^(l-r2)(l-r)  h3  =^(l-r)(l-5)-^(l-r2Xl-52) 

h2  =  ^(l-r)(l  +  s)--^(l-r2)(l-s2)  h4  =^(l  +  /-)(l-^)-^(l-r2)(l-^2) 

h5=(l-r2\l-s2) 


Figure  20.  Brick  Shape  Functions  Including  the  Bubble  Modes. 
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The  new  shape  function  will  increase  the  size  of  the  stiffness  matrix  to  a  ten  by 
ten  matrix.  This  means  that  there  are  now  ten  eigenvalues  and  ten  eigenvectors. 
The  original  eight  are  the  same,  but  the  two  new  ones  are  associated  with  the 
bubble  node.  These  are  called  bubble  modes  because  they  are  primarily  a 
displacement  of  the  bubble  node  without  much  deflection  in  the  other  modes. 

The  mode  shapes  remain  the  same  for  the  stiffness  matrix  including  the  bubble 
node.  The  only  change  beside  the  addition  of  the  bubble  modes  is  the  lowering 
of  the  bending  modes’  eigenvalue.  This  is  evidence  that  the  bending  stiffness  of 
this  element  is  reduced.  The  modes  and  their  eigenvalues  are  displayed  in  the 
figures  below. 


Shearing 

1=0.7575 


Uniaxial  Tension 
1=0.7575 


Biaxial  Tension 
1=1 .47 


Figure  21.  Brick  Element  Constant  Strain  Modes. 
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Figure  22.  Brick  Element  Bending  Modes. 
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Figure  23.  Brick  Element  Bubble  Modes. 
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Rotation  Horizontal  Translation  Vertical  Translation 

X=0  X=0  X=0 

Figure  24.  Brick  Element  Rigid  Body  Modes. 


6.2  Solution  Technique 

The  solution  starts  with  assuming  a  shape  function  of  the  solution.  The  solution 
is  quadratic  in  the  horizontal  and  the  vertical  coordinates,  as  well  as  their 
combination.  The  new  shape  functions  are  shown  below.  Notice  how  there  is  a 
quartic  polynomial  term.  It  is  included  in  the  same  manner  that  the  quadratic 
term  was  included  in  the  original  brick  development. 

u  =  Cj  +c2*  +  c3y  +  cAxy  +  cix2  +c6y2  +c7x2y  +  csxy2  +c9x2y2 

v  =  c10  +cnx  +  cny  +  cnxy  +  cux2  +c15y2  +c16.x:2y  +  c17xy2  +cKx2y1 

These  new  shape  functions  are  associated  with  eighteen  unknown  constants.  In 
the  same  manner  as  explained  in  the  section  on  why  the  brick  elements  lock,  the 
two  boundary  conditions  each  impose  three  constraints  on  each  of  the  new 
shape  functions.  Excluding  the  two  differential  equations,  there  are  twelve 
constraints.  This  leaves  six  constants  to  be  solved  for  by  the  differential 
equations.  The  equations  might  not  be  solved  for  exactly,  but  the  answer  can  be 
much  more  accurate  than  when  there  was  only  one  degree  of  freedom  that  could 
be  solved  for  in  the  two  differential  equations. 


Figure  25.  Bubble  Mode  Element  Deflection. 
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The  integration  of  the  stiffness  matrix  during  its  formulation  minimizes  the  error 
over  the  examined  structure.  This  produces  a  deformed  position  that  most 
accurately  enforces  all  the  boundary  conditions  as  well  as  the  differential 
equations.  The  linear  static  solution  with  the  mesh  described  using  fully 
integrated  bubble  modes  shows  a  tip  deflection  of  6.73865.  This  is  within  2.3% 
of  the  theoretical  solution,  a  vast  improvement  over  the  elements  without  the 
bubble  mode. 

6.3  Matrix  Inversion 

A  difficult  part  of  solving  a  linear  set  of  equations  is  that  there  can  be  n  unknowns 
in  n  equations,  where  n  can  be  very  large.  To  solve  for  any  particular  variable,  it 
is  necessary  to  have  one  equation  and  one  unknown.  Two  common  ways  of 
accomplishing  this  is  to  invert  the  matrix,  thereby  lumping  all  unknowns  into  one 
large  unknown,  and  substituting  equations  which  require  a  symbolic  processor, 
but  do  yield  a  series  of  one  unknown  equations. 

A  good  representation  of  the  second  part  of  the  substitution  method  is  shown  in 
the  following  equations.  Each  equation  only  uses  the  values  of  the  preceding 
equations  plus  an  additional  unknown.  The  first  two  equations  solved  for  the 
unknowns  are  shown  as  well.  The  pattern  is  evident  and  easily  repeatable. 
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Fmite  element  packages  use  this  process  to  solve  systems  of  equations  because 
they  are  much  faster  than  computing  the  inverse  of  a  large  matrix.  Unfortunately, 
they  must  solve  a  system  that  is  not  of  the  form  of  a  lower  triangular  matrix. 

There  is  a  way  to  represent  any  square  matrix  as  two  triangular  matrices  and  a 
diagonal  matrix,  and  a  symmetric  matrix  as  a  triangular  matrix  and  a  diagonal 
matrix.  This  method  is  shown  in  the  equations  below. 
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Now  a  single  fully  populated  linear  matrix  equation  can  be  solved  by  three  simple 
linear  equations. 


II 

SI 

LDLtu  =  F 

substituting 

y  =  DlJ  u 

Ly  =  F 

DlJ  u  =  y 

D.x  =  y 

substituting 

SI 

II 

Kl 

LT  u  =  x 

solve  for  unknown  vector 

This  might  seem  more  complicated  because  three  systems  of  equations  must  be 
solved  instead  of  a  single  system,  but  they  are  much  easier  systems  to  solve. 
They  also  have  an  added  benefit  when  dealing  with  a  highly  banded  matrix  as 
found  in  finite  element  applications.  None  of  the  matrices  lose  their  bandedness, 
unlike  the  inverse  of  a  banded  matrix.  The  inverse  of  a  banded  matrix  is 
normally  a  fully  populated  matrix  which  would  take  too  much  memory  to  perform 
calculations  on  it.  For  example,  a  ten  thousand-degree  of  freedom  system  would 
have  to  have  100  megabytes  of  memory  simply  to  hold  the  matrix  in  memory,  let 
alone  perform  operations  on  it.  The  simple  beam  model  shown  above  with  180 
elements  had  1260  degrees  of  freedom  and  would  have  needed  about  1 .6 
megabytes  of  memory  to  hold  the  inverted  matrix. 

6.4  Natural  Frequencies 

The  following  model  determined  the  natural  frequencies  of  the  system  using 
Patran.  It  is  the  same  180-element  model  described  above  comprised  of  Hex8 
brick  elements. 


Ea 

Figure  26.  Patran  Model  Used  to  Determine  Natural  Frequencies. 

The  following  table  shows  the  frequencies  determined  by  the  finite  element 
analysis.  The  accuracy  is  very  good  in  the  bending  modes  but  falls  off  in  the 
torsion  modes.  Partly  this  is  because  of  the  higher  strain  energy  in  those  modes. 
That  can  also  be  seen  in  the  higher  order  bending  modes,  although  those  have 
no  stress  gradient  in  the  x  direction. 
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FEM 

Calculations 

Error 

lat 

189.88 

194.52 

-2.39% 

tor 

94.28 

89.66 

5.16% 

tor 

288.27 

268.97 

7.18% 

tor 

497.92 

448.28 

1 1 .07% 

tor 

732.00 

627.60 

16.64% 

bend 

9.86 

9.73 

1.41% 

bend 

61.70 

60.95 

1 .22% 

bend 

173.13 

170.67 

1 .44% 

bend 

340.64 

334.44 

1 .85% 

bend 

565.60 

552.85 

2.31% 

The  mode  shapes  can  also  be  shown  as  part  of  the  results  in  finite  element 
analysis.  Below  are  several  of  the  mode  shapes  that  are  part  of  the  solution. 


Figure  27.  First  Bending  and  First  Lateral  Bending  Modes. 


Figure  28.  First  Torsion  and  Second  Bending  Modes. 


Figure  29.  Third  Bending  and  Second  Torsion  Modes. 

These  modes  compare  well  to  the  shapes  that  were  calculated  in  the  theorectical 
case.  The  enforcement  of  the  boundary  conditions,  both  the  natural  and 
geometric,  is  clear.  The  enforcement  of  the  zero  displacments  and  rotations  is 
visible  on  the  root  end.  The  zero  curvature  (zero  moment)  condition  can  be  seen 
in  each  mode  on  the  free  end.  It  is  harder  to  see  the  zero  shear  condition,  but 


6-6 


the  free  face  of  the  beam  remains  perpendicular  to  the  edges  leading  to  the  end 
face. 

6.5  Element  Applicability  and  Meshing  Criteria 

It  is  interesting  to  note  that  the  model  run  with  the  shell  elements  with  reduced 
order  integration  had  more  accurate  results  with  fewer  elements.  This  is 
because  the  theory  of  plates  is  based  on  the  same  stress  strain  relationship  that 
the  shell  elements  use. 

Solid  elements  are  much  more  useful  when  there  is  a  nonuniform  thickness  or  if 
there  is  a  dramatic  stress  gradient  through  the  thickness  caused  by  forces  other 
than  bending.  It  is  much  easier  to  correctly  vary  the  thickness  of  the  solid 
elements  along  a  span  than  on  the  shell  elements.  Additionally,  it  is  possible  to 
put  more  than  a  single  layer  of  elements  through  the  thickness  unlike  in  shell 
elements. 

Although  the  element  is  no  longer  as  susceptible  to  being  overly  stiff,  this  does 
not  mean  that  the  aspect  ratio  of  the  elements  does  not  matter  at  all.  It  is  still 
important  to  keep  the  brick  as  close  to  a  perfect  cube  as  practical  for  a  given 
situation.  This  is  because  multiplying  the  formulation  by  the  determinate  of  the 
Jacobian  forms  the  stiffness  matrix.  If  the  element  is  very  skewed,  this 
transformation  into  the  local  coordinates  is  not  as  accurate.  To  keep  the  error  to 
a  minimum  the  skew  of  element  must  be  at  a  minimum  practical  level.  In  our 
example  it  would  be  best  to  have  all  elements  a  cube  0.1  on  each  edge.  This 
would  unnecessarily  increase  the  mesh  to  2000  elements.  As  long  as  the 
element  edges  remain  normal  the  aspect  ratio  can  be  left  higher  than  1 . 
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SECTION  7 

BRICK  ELEMENT  WITH  REDUCED  INTEGRATION 

The  final  element  that  we  will  look  at  is  the  brick  element  with  reduced 
integration.  There  is  no  bubble  mode,  so  the  element  is  comparable  to  the 
reduced  integration  shell  element.  The  element  still  uses  the  brick  element 
definition  of  the  strain  energy. 

The  same  formulations  of  the  strain  energy  and  stress-strain  relationship  are 
used.  Just  as  we  used  a  single  integration  point  in  the  reduced  shell  integration, 
there  will  only  be  one  integration  point  in  each  of  the  local  coordinate  directions. 

Mode  shapes  are  not  changed  by  reduced  integration,  although  the  associated 
eigenvalues  are.  These  new  eigenvalues  dramatically  change  the  energy 
required  to  deform  the  structure,  bringing  the  solution  towards  the  exact  solution. 

7.1  Mode  Shapes 

The  reduced  order  integration  technique  changes  the  associated  eigenvalues  for 
some  mode  shapes.  These  are  caused  by  the  modification  of  the  stiffness  matrix 
when  reduced  order  Gauss  point  integration  incorrectly  evaluates  it.  The  theory 
is  the  same  as  it  is  in  the  shell  elements;  the  changes  allow  the  correct  mode  to 
deform  with  less  energy  being  used  by  the  same  deflection. 

This  theory  works,  but  not  as  well  as  in  the  shell  element.  The  main  reason  that 
it  is  worse  than  the  shell  element  is  that  the  eigenvalue  is  reduced  to  zero  instead 
of  some  new  positive  value.  Both  were  similar  in  that  only  the  bending  modes 
were  effected.  The  goal  was  to  reduce  the  energy  required  to  excite  these 
modes,  so  the  method  worked  well  in  that  respect. 

7.2  Deflection  Solution 

The  change  is  apparent  in  the  solution.  The  solution  quickly  approaches  the 
closed  form  solution.  An  interesting  thing  happens  as  the  number  of  elements  is 
reduced  to  just  a  couple  elements.  The  zero  eigenvalues  of  the  bending  modes 
prevent  the  stiffness  matrix  from  being  inverted.  After  a  few  of  the  elements  are 
connected  the  complete  system  can  no  longer  exhibit  a  single  bending  mode,  so 
the  matrix  is  no  longer  singular.  The  individual  elements  still  have  zero  energy 
associated  with  their  bending  so  the  solution  will  not  be  exact.  All  strain  energy 
is  based  on  the  other  elastic  modes  although  the  solution  shows  primarily 
bending.  The  convergence  is  close  to  the  exact  solution,  but  it  is  not  as  close  as 
the  bubble  mode  solution. 

If  the  elements  do  not  converge  to  the  correct  solution,  they  do  converge  quickly. 
The  difference  between  30  spanwise  elements  (180  total  elements)  and  50 
spanwise  elements  (500  total)  is  negligible. 


7-1 


Figure  30.  Reduced  Integration  Defection  Solution 


7.3  Natural  Frequencies 

There  is  also  an  associated  system  frequency  change  with  the  new  elements. 
Much  like  the  deflection  solution,  the  frequencies  come  closer  to  the  exact 
solution,  but  not  as  close  as  the  bubble  functions.  Depending  on  the  application 
the  frequencies  may  not  be  close  enough  to  accept.  They  still  have  about  ten- 
percent  error  in  the  low  frequency  modes. 


FEM 

Calculations 

Error 

1 

lat 

190.43 

194.52 

-2.10% 

1 

tor 

94.88 

89.66 

5.83% 

2 

tor 

291.17 

268.97 

8.25% 

3 

tor 

506.14 

448.28 

12.91% 

4 

tor 

749.94 

627.60 

19.49% 

1 

bend 

10.59 

9.73 

8.88% 

2 

bend 

66.16 

60.95 

8.54% 

3 

bend 

186.05 

170.67 

9.01% 

4 

bend 

367.41 

334.44 

9.86% 

5 

bend 

612.33 

552.85 

10.76% 
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The  error  in  frequency  stems  from  the  same  cause  as  the  error  in  the  deflection. 
The  system  has  the  incorrect  energy  associated  with  the  given  deflection  modes, 
changing  the  amount  of  energy  (frequency)  required  to  excite  a  given  mode. 

7.4  Comparison  with  Bubble  Node 

The  solution  is  not  as  accurate  as  the  bubble  mode  in  this  solution.  That  is 
because  instead  of  making  the  model  more  accurate  by  increasing  the  order  of 
the  shape  functions  to  better  enforce  the  strain  displacement  relationship,  the 
most  prominent  deflection  mode  has  the  eigenvalue  arbitrarily  reduced.  It  does 
make  sense  that  if  a  system  should  bend  and  it  doesn’t,  reducing  energy 
associated  with  the  bending  mode  will  make  the  system  approach  the  correct 
value.  The  problem  is  that  a  better  solution  comes  from  the  degradation  of  the 
model  and  not  from  a  more  exact  model. 

Another  way  to  look  at  the  solution  is  to  say  that  the  shape  functions  are  a  series 
of  normal  functions  added  together.  As  the  number  of  normal  functions 
increases,  the  solution  approaches  the  exact  solution  as  in  the  Fourier  Series 
expansion  of  a  function.  It  can  be  shown  that  increasing  the  polynomial  order  of 
the  shape  function  will  bring  the  solution  closer  to  the  exact  solution.  That  is  the 
theory  behind  p-elements,  which  increase  the  order  of  the  polynomial  shape 
function  until  the  solution  converges  within  a  given  tolerance. 

Conversely,  the  reduced  order  integration  has  no  proof  that  the  mesh  will  even 
converge  to  the  right  solution  let  alone  a  more  accurate  one  than  a  fully 
integrated  one.  In  practice  it  does  converge  towards  the  correct  solution  most  of 
the  time,  but  the  analyst  should  be  aware  that  there  is  no  guarantee  that  it  will. 

The  reason  that  the  bending  mode  eigenvalue  goes  to  zero  using  reduced 
integration  is  that  the  single  Gauss  point  is  at  a  zero  strain  position  in  the 
element.  That  means  that  it  will  have  zero  energy  and  a  zero  eigenvalue.  This 
can  be  seen  in  the  figure  below. 

The  single  Gauss  point  in  the 
center  of  the  element  is 
located  in  the  middle  of 
symmetric  deflection  so 
there  is  no  strain  at  that 
point. 


Figure  31.  Gauss  Point  Location  in  Reduced  Order  Bending  Mode. 

In  the  figure  you  can  also  see  the  locations  of  the  fully  integrated  Gauss  points  by 
the  intersection  of  the  dashed  lines.  It  is  clear  that  they  are  on  lines  that  deform, 
and  therefore  have  an  energy  associate  with  that  deflection. 
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SECTION  8 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  goal  of  any  structural  analysis  is  to  correctly  predict  the  conditions  that  will 
exist  in  a  given  structure.  The  best  way  to  accomplish  this  when  using  these 
elements  that  are  subject  to  being  over  stiff  is  to  understand  what  the  different 
element  property  choices  mean  mathematically  in  the  model. 

The  reduced  integration  technique  works  well  in  the  shell  method  because  the 
bending  strain  definition  is: 

El(ff 

and  phi  prime  is  a  non-zero  constant  at  the  Gauss  point.  This  implies  that  there 
will  be  a  non-zero  eigenvalue  associated  with  the  bending  mode  because  there  is 
strain  at  the  Gauss  point. 

In  the  brick  element,  the  Gauss  point  is  at  a  saddle  point.  The  values  of  the 
gradients  in  both  the  r  and  the  s  directions  are  zero  at  the  Gauss  point,  as  are  the 
shear  gradients.  This  means  that  the  strain  at  the  Gauss  point  is  zero  and  will 
produce  a  zero  eigenvalue.  The  saddle  point  can  be  seen  in  the  figures  showing 
u  deflection  in  a  bending  mode. 


Figure  32.  Deflection  Plot  of  Bending  Mode  Shape. 


The  main  difference  is  again  due  to  the  difference  in  the  strain  energy  definition. 
The  Cauchy  formulation  used  in  the  brick  element  is  not  accurate  in  these 
conditions  dominated  by  bending  and  shear. 

The  reduced  order  shell  element  is  the  best  solution  for  overcoming  locking  in  the 
shell  element  because  the  stiffness  of  the  bending  mode  is  not  reduced  to  zero. 
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It  is  the  best  formulation  of  the  element  to  use  in  most  conditions.  The  analyst 
should  still  be  aware  that  there  is  still  no  method  to  know  that  the  solution  that  the 
model  converges  to  is  the  exact  solution,  but  in  most  cases  it  is  closer  to  reality 
than  the  fully  integrated  element. 

The  best  brick  element  formulation  to  use  is  the  element  that  has  quadratic 
shape  functions.  That  means  that  the  simplest  element  would  be  the  eight-node 
element  using  bubble  functions.  Other  elements  that  also  use  the  quadratic 
shape  functions  are  the  20,  21 , 26,  and  27  node  elements.  These  will 
dramatically  increase  the  computational  time  if  the  number  of  elements  remains 
constant. 

Modern  FEA  code  provides  hourglass  mode  stiffening.  This  is  a  change  in  the 
stiffness  matrix  that  can  provide  stiffness  to  those  modes  that  normally  would 
have  a  zero  eigenvalue  under  reduced  integration.  This  is  another  solution  to  the 
problems  associated  with  reduced  integration  of  brick  elements. 

Whatever  formulation  and  elements  used  in  the  FEM,  it  is  very  important  that  the 
person  using  them  knows  how  and  why  they  act  as  they  do.  An  analysis  is  only 
as  good  as  the  method  used,  so  knowing  that  you  are  using  the  right  method  for 
the  particular  analysis  is  extremely  important. 
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